#ifndef __FIRSTORDER_H__
#define __FIRSTORDER_H__

#include "reduce.h"
#include "framework.h"
#include "kernel.h"
#include <float.h>

__device__ void coopExtractRowVector(float* data, int dataPitch, int index, int dimension, float* destination) {
  float* xiRowPtr = data + (index * dataPitch) + threadIdx.x;
  for(int currentDim = threadIdx.x; currentDim < dimension; currentDim += blockDim.x) {
    destination[currentDim] = *xiRowPtr;
    xiRowPtr += blockDim.x;
  }
}

int firstOrderPhaseOneSize(bool iLowCompute, bool iHighCompute, int n_dim) {
	int size = 0;
	if (iHighCompute) { size+=sizeof(float)*n_dim; }
	if (iLowCompute)  { size+=sizeof(float)*n_dim; }
	return size;
}
int firstOrderPhaseTwoSize() {
	int size = 0;
	return size;
}

template<bool iLowCompute, bool iHighCompute>
  __global__ void	firstOrderPhaseOne(float* devData, int devDataPitchInFloats, float* devTransposedData, int devTransposedDataPitchInFloats, float* devLabels, int n_obs, int n_dim, float epsilon, float cost_epsilon, float* devAlpha, float* devF, float alpha1Diff, float alpha2Diff, int iLow, int iHigh, float* devCache, int devCachePitchInFloats, int iLowCacheIndex, int iHighCacheIndex, int* devLocalIndicesRL, int* devLocalIndicesRH, float* devLocalFsRL, float* devLocalFsRH) {

    //
    extern __shared__ float xIHigh[];
	float* xILow;
	__shared__ int tempLocalIndices[BLOCKSIZE];
	__shared__ float tempLocalFs[BLOCKSIZE];
  
  
    if (iHighCompute) {
		xILow = &xIHigh[n_dim];
	} else {
		xILow = xIHigh;
	}
  
    if (iHighCompute) {
        //Load xIHigh into shared memory
        coopExtractRowVector(devTransposedData, devTransposedDataPitchInFloats, iHigh, n_dim, xIHigh);
	}
    if (iLowCompute) {
        //Load xILow into shared memory
        coopExtractRowVector(devTransposedData, devTransposedDataPitchInFloats, iLow, n_dim, xILow);
    }
	__syncthreads();
	
	int globalIndex = IMUL(blockDim.x, blockIdx.x) + threadIdx.x;

    float alpha;
    float f;
    float label;
    int reduceFlag;
  
	if (globalIndex < n_obs) {
        alpha = devAlpha[globalIndex];
        f = devF[globalIndex];
    	label = devLabels[globalIndex];
     
		if (alpha > epsilon) {
			if (alpha < cost_epsilon) {
				reduceFlag = REDUCE0 | REDUCE1; //Unbound support vector (I0)
			} else {
				if (label > 0) {
					reduceFlag = REDUCE0; //Bound positive support vector (I3)
				} else {
					reduceFlag = REDUCE1; //Bound negative support vector (I2)
			}}
		} else {
			if (label > 0) {
				reduceFlag = REDUCE1; //Positive nonsupport vector (I1)
			} else {
				reduceFlag = REDUCE0; //Negative nonsupport vector (I4)
		}}
	} else {
		reduceFlag = NOREDUCE;
	}

  
    float highKernel = 0;
    float lowKernel = 0;
	if (reduceFlag > 0) {
		if (!iHighCompute) {
			highKernel = devCache[(devCachePitchInFloats * iHighCacheIndex) + globalIndex];
		}
		if (!iLowCompute) {
			lowKernel = devCache[(devCachePitchInFloats * iLowCacheIndex) + globalIndex];
		}
        if (iHighCompute && iLowCompute) {
            dualKernel(devData+globalIndex,devDataPitchInFloats,devData+globalIndex+(devDataPitchInFloats*n_dim),xIHigh,1,xILow,1,highKernel,lowKernel);
        } else if (iHighCompute) {
            highKernel = kernel(devData + globalIndex, devDataPitchInFloats, devData + globalIndex + (devDataPitchInFloats * n_dim), xIHigh, 1);
        } else if (iLowCompute) {
            lowKernel = kernel(devData + globalIndex, devDataPitchInFloats, devData + globalIndex + (devDataPitchInFloats * n_dim), xILow, 1);
        }

		f = f + alpha1Diff * highKernel;
		f = f + alpha2Diff * lowKernel;

        if (iLowCompute) {
			devCache[(devCachePitchInFloats * iLowCacheIndex) + globalIndex] = lowKernel;
		}
		if (iHighCompute) {
			devCache[(devCachePitchInFloats * iHighCacheIndex) + globalIndex] = highKernel;
		}

		devF[globalIndex] = f;
	}
	__syncthreads();

    
    // BEGIN REDUCTION PHASE
    if ((reduceFlag & REDUCE0) == 0) {
        tempLocalFs[threadIdx.x] = -FLT_MAX; //Ignore me
    } else {
        tempLocalFs[threadIdx.x] = f;
        tempLocalIndices[threadIdx.x] = globalIndex;
    }
    __syncthreads();
  
    argmaxReduce(tempLocalFs, tempLocalIndices);
    
    if (threadIdx.x == 0) {
        devLocalIndicesRL[blockIdx.x] = tempLocalIndices[0];
		devLocalFsRL[blockIdx.x] = tempLocalFs[0];
    }
    __syncthreads();

    if ((reduceFlag & REDUCE1) == 0) {
        tempLocalFs[threadIdx.x] = FLT_MAX; //Ignore me
    } else {
        tempLocalFs[threadIdx.x] = f;
        tempLocalIndices[threadIdx.x] = globalIndex;
    }
    __syncthreads();
   
    argminReduce(tempLocalFs, tempLocalIndices);
    if (threadIdx.x == 0) {
		devLocalIndicesRH[blockIdx.x] = tempLocalIndices[0];
		devLocalFsRH[blockIdx.x] = tempLocalFs[0];
    }

    //end first phase
}


__global__ void firstOrderPhaseTwo(float* devData, int devDataPitchInFloats, float* devTransposedData, int devTransposedDataPitchInFloats, float* devLabels, float* devKernelDiag, float* devAlpha, void* devResult, float cost, int n_dim, int* devLocalIndicesRL, int* devLocalIndicesRH, float* devLocalFsRL, float* devLocalFsRH, int inputSize) {
  __shared__ int tempIndices[BLOCKSIZE];
  __shared__ float tempFs[BLOCKSIZE];
  
  
  //Load elements
  if (threadIdx.x < inputSize) {
    tempIndices[threadIdx.x] = devLocalIndicesRH[threadIdx.x];
    tempFs[threadIdx.x] = devLocalFsRH[threadIdx.x];
  } else {
    tempFs[threadIdx.x] = FLT_MAX;
  }

  if (inputSize > BLOCKSIZE) {
    for (int i = threadIdx.x + BLOCKSIZE; i < inputSize; i += blockDim.x) {
      argMin(tempIndices[threadIdx.x], tempFs[threadIdx.x], devLocalIndicesRH[i], devLocalFsRH[i], tempIndices + threadIdx.x, tempFs + threadIdx.x);
    }
  }
  __syncthreads();
  argminReduce(tempFs, tempIndices);
  int iHigh = tempIndices[0];
  float bHigh = tempFs[0];

  //Load elements
  if (threadIdx.x < inputSize) {
    tempIndices[threadIdx.x] = devLocalIndicesRL[threadIdx.x];
    tempFs[threadIdx.x] = devLocalFsRL[threadIdx.x];
  } else {
    tempFs[threadIdx.x] = -FLT_MAX;
  }

  if (inputSize > BLOCKSIZE) {
    for (int i = threadIdx.x + BLOCKSIZE; i < inputSize; i += blockDim.x) {
      argMax(tempIndices[threadIdx.x], tempFs[threadIdx.x], devLocalIndicesRL[i], devLocalFsRL[i], tempIndices + threadIdx.x, tempFs + threadIdx.x);
    }
  }
  __syncthreads();
  
  argmaxReduce(tempFs, tempIndices);
  
  int iLow = tempIndices[0];
  float bLow = tempFs[0];

  
  
  float* highPointer = devTransposedData + (iHigh * devTransposedDataPitchInFloats);
  float* lowPointer = devTransposedData + (iLow * devTransposedDataPitchInFloats);  
  
  parallelKernel(highPointer, highPointer + n_dim, lowPointer, tempFs);

  if (threadIdx.x == 0) {
    
    float eta = devKernelDiag[iHigh] + devKernelDiag[iLow];
      
    float kernelEval = tempFs[0];
    eta = eta - 2*kernelEval;
      
    float alpha1Old = devAlpha[iHigh];
    float alpha2Old = devAlpha[iLow];
    float alphaDiff = alpha2Old - alpha1Old;
    float lowLabel = devLabels[iLow];
    float sign = devLabels[iHigh] * lowLabel;
    float alpha2UpperBound;
    float alpha2LowerBound;
    if (sign < 0) {
      if (alphaDiff < 0) {
        alpha2LowerBound = 0;
				alpha2UpperBound = cost + alphaDiff;
      } else {
        alpha2LowerBound = alphaDiff;
        alpha2UpperBound = cost;
      }
    } else {
      float alphaSum = alpha2Old + alpha1Old;
      if (alphaSum < cost) {
        alpha2UpperBound = alphaSum;
        alpha2LowerBound = 0;
      } else {
        alpha2LowerBound = alphaSum - cost;
        alpha2UpperBound = cost;
			}
    }
    float alpha2New;
    if (eta > 0) {
      alpha2New = alpha2Old + lowLabel*(bHigh - bLow)/eta;
      if (alpha2New < alpha2LowerBound) {
        alpha2New = alpha2LowerBound;
      } else if (alpha2New > alpha2UpperBound) {
				alpha2New = alpha2UpperBound;
      }
    } else {
      float slope = lowLabel * (bHigh - bLow);
      float delta = slope * (alpha2UpperBound - alpha2LowerBound);
      if (delta > 0) {
        if (slope > 0) {
          alpha2New = alpha2UpperBound;
        } else {
          alpha2New = alpha2LowerBound;
        }
      } else {
        alpha2New = alpha2Old;
      }
    }
    float alpha2Diff = alpha2New - alpha2Old;
    float alpha1Diff = -sign*alpha2Diff;
    float alpha1New = alpha1Old + alpha1Diff;

    *((float*)devResult + 0) = alpha2Old;
    *((float*)devResult + 1) = alpha1Old;
    *((float*)devResult + 2) = bLow;
    *((float*)devResult + 3) = bHigh;
    devAlpha[iLow] = alpha2New;
    devAlpha[iHigh] = alpha1New;
    *((float*)devResult + 4) = alpha2New;
    *((float*)devResult + 5) = alpha1New;
    *((int*)devResult + 6) = iLow;
    *((int*)devResult + 7) = iHigh;
  }
}


void launchFirstOrder(bool iLowCompute, bool iHighCompute, int n_obs, int n_dim, dim3 blocks, dim3 threads, dim3 globalThreads, float* devData, int devDataPitchInFloats, float* devTransposedData, int devTransposedDataPitchInFloats, float* devLabels, float epsilon, float cost_epsilon, float* devAlpha, float* devF, float sAlpha1Diff, float sAlpha2Diff, int iLow, int iHigh, float* devCache, int devCachePitchInFloats, int iLowCacheIndex, int iHighCacheIndex, int* devLocalIndicesRL, int* devLocalIndicesRH, float* devLocalFsRH, float* devLocalFsRL, float* devKernelDiag, void* devResult, float cost) {

  int phaseOneSize = firstOrderPhaseOneSize(iLowCompute, iHighCompute, n_dim);
  int phaseTwoSize = firstOrderPhaseTwoSize();
  if (iLowCompute == true) {
    if (iHighCompute == true) {
        firstOrderPhaseOne<true, true><<<blocks, threads, phaseOneSize>>>(devData, 
                devDataPitchInFloats, devTransposedData, devTransposedDataPitchInFloats, devLabels, n_obs, 
                n_dim, epsilon, cost_epsilon, devAlpha, devF, sAlpha1Diff, sAlpha2Diff, iLow, iHigh, 
                devCache, devCachePitchInFloats, iLowCacheIndex, iHighCacheIndex, devLocalIndicesRL, 
                devLocalIndicesRH, devLocalFsRL, devLocalFsRH);
    } else if (iHighCompute == false) {
       firstOrderPhaseOne<true, false><<<blocks, threads, phaseOneSize>>>(devData, 
               devDataPitchInFloats, devTransposedData, devTransposedDataPitchInFloats, devLabels, n_obs, 
               n_dim, epsilon, cost_epsilon, devAlpha, devF, sAlpha1Diff, sAlpha2Diff, iLow, iHigh, 
               devCache, devCachePitchInFloats, iLowCacheIndex,iHighCacheIndex, devLocalIndicesRL, 
               devLocalIndicesRH, devLocalFsRL, devLocalFsRH);
    } 
  } else if (iLowCompute == false) {
    if (iHighCompute == true) {
       firstOrderPhaseOne<false, true><<<blocks, threads, phaseOneSize>>>(devData, 
               devDataPitchInFloats, devTransposedData, devTransposedDataPitchInFloats, devLabels, n_obs, 
               n_dim, epsilon, cost_epsilon, devAlpha, devF, sAlpha1Diff, sAlpha2Diff, iLow, iHigh, devCache, 
               devCachePitchInFloats, iLowCacheIndex, iHighCacheIndex, devLocalIndicesRL, devLocalIndicesRH, 
               devLocalFsRL, devLocalFsRH);
    } else if (iHighCompute == false) {
       firstOrderPhaseOne<false, false><<<blocks, threads, phaseOneSize>>>(devData, 
               devDataPitchInFloats, devTransposedData, devTransposedDataPitchInFloats, devLabels, n_obs, 
               n_dim,epsilon,cost_epsilon, devAlpha, devF, sAlpha1Diff, sAlpha2Diff, iLow, iHigh, devCache, 
               devCachePitchInFloats,iLowCacheIndex, iHighCacheIndex, devLocalIndicesRL, devLocalIndicesRH, 
               devLocalFsRL, devLocalFsRH);
    }
  }
 
  //run phase 2 
  firstOrderPhaseTwo<<<1, globalThreads, phaseTwoSize>>>(devData, devDataPitchInFloats, 
          devTransposedData, devTransposedDataPitchInFloats, devLabels, devKernelDiag, devAlpha, devResult, 
          cost, n_dim, devLocalIndicesRL, devLocalIndicesRH, devLocalFsRL, 
          devLocalFsRH, blocks.x);
}
#endif
